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We report calculation of heat capacity of an attractive Bose-Einstein condensate, with the number 
N of bosons increasing and eventually approaching the critical number N cr for collapse, using the 
correlated potential harmonics (CPH) method. Boson pairs interact via the realistic van der Waals 
potential. It is found that the transition temperature T c increases initially slowly, then rapidly as N 
becomes closer to N cr . The peak value of heat capacity for a fixed N increases slowly with N, for 
N far away from N cr . But after reaching a maximum, it starts decreasing when N approaches N cr . 
The effective potential calculated by CPH method provides an insight into this strange behavior. 
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I. INTRODUCTION 

Bose-Einstein condensation (BEC) is the transition 
process, in which a macroscopic fraction of bosons 
goes into the lowest energy state, as the temperature 
is lowered below a certain critical temperature T c [l|. 
It was predicted by Einstein in 1925, based on Bosc's 
explanation of black body radiation. A great deal of 
activity, both theoretical and experimental, has been 
seen in this field, since the experimental realization of 
BEC in 1995. Although a number of static, dynamic 
and thermodynamic properties have been studied 0, [3] , 
not much attention has been paid to the heat capacity 
of attractive condensates. The main motivation of this 
work is to fill this gap. 

In laboratory experiments, the condensate is trapped 
by a confining potential, usually a harmonic oscillator 
potential. An attractive condensate (e.g. 7 Li con- 
densate) has a negative value of the s-wave scattering 
length a s and collapses, when the number of particles 
N in the condensate exceeds a critical number N cr . 
On the other hand a repulsive condensate (e.g. 87 Rb 
condensate) corresponds to a s > and is stable for any 
N, since repulsively interacting bosons are contained 
in the externally applied trap. The situation is quite 
different for an attractive condensate: attractive bosons 
tend to come to the center of the trap, which is balanced 
only by the kinetic pressure, resulting in a metastablc 
condensate. The total attraction increases as the number 
of pairs N(N — l)/2, while the kinetic pressure increases 
as N. Thus for N larger than a critical value N cri the 
net attraction dominates and a collapse occurs. 

In this communication, we report the calculation of 
heat capacity of an attractive condensate containing a 
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fixed number of 7 Li atoms, interacting via the realistic 
van der Waals potential, appropriate for the experi- 
mental scattering length. The features are markedly 
different from those of repulsive condensates, only which 
have so far been investigated. In a repulsive condensate, 
the heat capacity Cjv(T) for a fixed number N of bosons 
in the trap, as also the critical temperature T c ^moothly 
approach a constant value as N — s- oo 0, Q. As a 
function of T, the heat capacity for a given N increases 
to a maximum (Cjv)toox; then falls rapidly to a satura- 
tion value 3NkB as T increases (fee is the Boltzmann 
constant). These features are qualitatively similar to 
those of a trapped non- interacting condensate [!, [f|. 
However for an attractive condensate, there are impor- 
tant changes in the nature. This is due to the fact that 
the number of available energy levels of the system is 
limited, especially when N ~ > N cr , while for any N, 
there are infinitely many energy levels for the repulsive 
or non-interacting condensate. In the limit of high T, 
both the repulsive and the non-interacting condensate 
behave as the corresponding trapped Bosc-gas, resulting 
in a saturation in Cjy(T). On the other hand, an 
attractive condensate also shows similar behavior, only 
if it is allowed to absorb energy internally through ro- 
tational motion involving large orbital angular momenta. 



We provide an understanding of this peculiar nature 
based on the many-body picture. For the theoretical 
calculation, we adopt the correlated potential harmonic 
(CPH) method [f| Q to solve the many-body problem 
approximately. This technique is based on the potential 
harmonics (PH) expansion method @. The laboratory 
BEC must be very dilute to preclude three-body colli- 
sions, which lead to molecule formation and consequent 
depletion. Hence only two-body correlations are rele- 
vant. The PH is a subset [8| of the full hyperspherical 
harmonics (HH) basis @, that involves only two-body 
correlations. Hence the PH basis is a good approxima- 
tion for expanding the condensate wave function. It 
reduces the bulk of the numerical procedure immensely, 
while retaining the most important basic features of the 
condensate. However, the leading members of the PH 
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basis do not have the correct short separation behavior 
of the interacting Faddeev component. This causes 
a very slow rate of convergence of the PH expansion 
basis. To correct for this, we include a short-range 
correlation function in the expansion basis. This corre- 
lation function is obtained as the zero-energy solution 
of the two-body Schrodinger equation Q . It is a correct 
representation of the short separation behavior and also 
incorporates the s-wave scattering length, a s . through 
its asymptotic behavior @. The technique has been 
shown to reproduce known results, both experimental 
and theoretical ■ These include the following: ground 
state properties (energy, wave function, condensate 
size, one-body density, pair-correlation, etc.) as also 
multipolar moments of both repulsive and attractive 
condensates, correct prediction of the critical number 
and collapse scenario of attractive condensates, thermo- 
dynamic properties of repulsive condensates, properties 
of condensates in finite traps, etc. 

We can understand the behavior of heat capacity 
of attractive condensates in terms of the energy levels 
of the system produced by the CPH method. This 
method generates an effective potential in which the 
condensate moves. For an attractive condensate, the 
effective potential has a metastable region (MSR), 
separated from a deep well on the inner side by an 
intermediate finite barrier. A finite number of energy 
levels are supported by the MSR. As the number N 
of atoms increases, the MSR shrinks and the number 
of energy levels reduce drastically. As temperature in- 
creases, particles are distributed in higher energy levels, 
according to Bose distribution formula. Thus at low 
temperatures the internal energy and CV(T) increase 
with temperature. At higher temperatures, the bosons 
have fewer levels to occupy, causing Cn {T) to differ from 
the repulsive case. There is also a dominant effect of the 
drastically reducing number of energy levels as N — > N cr . 

The paper is organized as follows. For easy readability 
and to introduce our notations, we briefly review the 
correlated potential harmonic method in Section II. 
Section III provides our numerical procedure. Results 
and discussion arc presented in Section IV. Finally we 
draw our conclusions in Section V. 
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where Xi is the position vector of the i-th particle. The 
Schrodingcr equation governing the relative motion of the 
system trapped in a harmonic well, is 



— £v| +F trap (Ci,...,Cv) 



V(Ci,.~,bf)-E R \l)(Ci,...,bf) = 0, (2) 

where j\f = N — 1 and the trapping potential V trap and 
the interatomic interaction V are expressed in terms of 
the Jacobi vectors. The energy of the relative motion is 
En. Next, we introduce hyperspherical variables corre- 
sponding to the set of j\f Jacobi vectors. First, a hypcr- 
radius is defined as 



.2 = 1 



(3) 



The remaining set of (3J\f — 1) 'hyperangles' consists of 
27V polar angles of TV Jacobi vectors and (M — 1) angles 
defining their relative lengths @. In the hyperspherical 
harmonics expansion method (HHEM) ip is expanded 
in the complete set of hyperspherical harmonics (HH), 
which are the eigenfunctions of the grand orbital oper- 
ator [hypcrangular part of the j\f dimensional Laplace 
operator, given by the sum in the first term of Eq. 
(0)] 0. Substitution of this in Eq. ^ and projection 
on a particular HH result in a set of coupled differential 
equations. Imposition of symmetry of the wave function 
and calculation of the matrix elements become increas- 
ingly difficult and tedious as N increases. In addition, 
the degeneracy of the HH basis increases very rapidly @ 
with the increase in the grand orbital quantum number 
K, Hence a convergent calculation using HHEM with 
the full HH basis is extremely computer intensive and 
unmanageable for N > 3. This is the price one pays for 
keeping all many-body correlations in tp. 



II. CORRELATED POTENTIAL HARMONICS 
(CPH) METHOD 

We adopt the correlated potential harmonics 
method @, Q to solve the many-body problem of 
the BEC. We briefly recapitulate the technique in the 
following. Interested readers can find details in the cited 
references. 

For the relative motion of a system of N identical spin- 



However all these complications can be avoided and 
a much simpler computational procedure can be formu- 
lated for the laboratory BEC, which is designed to be 
extremely dilute (typical number density is ~ I0 15 cm~ 3 ) 
in order to avoid recombination through three-body col- 
lisions. Thus three-body correlations and three-body 
forces are totally negligible. We can then express ip as 
a sum of two-body Faddeev component ipij for the (ij)- 
interacting pair 
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Note that the assuption of two-body correlations only 
makes ipij a function of the pair separation vector and 
the hypcrradius only. One can then expand ipij in a sub- 
set of HH, called the potential harmonics (PH) subset, 
which is sufficient for the expansion of the interaction 
potential V"(fy) as a function in the hyperangular space 
for the (repartition. Since the labeling of the particles 
is arbitrary, we can choose = £w- Then the corre- 
sponding PH, V^K+ii^Af) (t ne argument is the full set 
of hypcranglcs for the (repartition) is independent of 

{Ci; • ■ • 7 C/V-i} and a simple analytic expression is possi- 
ble @. Expansion of the Faddeev component in the PH 
basis reads 



ipij{r i: 
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where ^ff+z^jv") is a potential harmonic The in- 
dependent factor in front is included to remove the first 
derivative with respect to r. Substitution of this expan- 
sion in the Faddeev equation for the (ij)-partition 
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[where T = -^Ei=i^-] an< ^ projection on the PH 
corresponding to the (repartition give a set of cou- 
pled differential equations in r. Note that any realistic 
two-body potential, V(fij) can be used. A realistic in- 
teratomic potential has a very strong repulsion (arising 
from the nucleus- nucleus repulsion) at very short separa- 
tions. Consequently, corresponding ipij must be vanish- 
ingly small for small values of nj. But the leading PH 
(corresponding to K = 0) in the expansion in Eq. ([5]) is 
a constant and does not have this behavior. Hence con- 
vergence of the expansion in Eq. ([5]) will be very slow. 
To improve the rate of convergence, we include a short- 
range correlation function rj{rij) in the expansion basis, 
so that Eq. is replaced by 



ipij(rij,r) = r 
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The short-range correlation function is chosen to have the 
expected behavior of (/', ; (r) y , r) for small r- n in the fol- 
lowing manner. The small ry behavior of ipij will be that 
of a zero-energy pair interacting via V(fij), since the en- 
ergy of the interacting pair is practically zero. We obtain 
r](fij) by solving the zero-energy two-body Schrodinger 
equation 



h 2 1 d 

m rjj dnj 



dm 



+ V(nj)r,(nj) = o. (8) 



Inclusion of the short-range correlation function rj(fij) 
enhances the rate of convergence greatly, which has been 
checked in our numerical calculation. 



The laboratory BEC is very dilute; hence the average 
separation of the atoms is very large compared with 
the range of interatomic interactions. Moreover, the 
atoms scatter with almost zero energy. Hence the 
effective two-body interaction is represented by the 
s-wave scattering length a s . In our calculation, we take 
V(fij) to be the van der Waals potential with a hard 
core: V(fij) = — -f^ for r,-j > r c and = oo for < r c . 
The correlation function obtained by solving Eq. © 
quickly attains its asymptotic form C(l — f^-) for large 
Tij . The asymptotic normalization is chosen to make the 
wavefunction positive at large . The hard core radius 
r c is adjusted so that the calculated a s is the actual 
experimental value of the scattering length 0). This 
procedure assures that the realistic two-body interaction 
appropriate for the condensate has been incorporated. 

Substitution of the expansion, Eq. (JT]) in Eq. (O and 
projection on the PH corresponding to the (r/)-partition 
result in 
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+ Vtra P (r) - E R \U K i(r) 

+ ^2fKiVKK'{r)fK'iU K n(r) = 0, 



(9) 
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where U K i{r) = fi<iu l K (r), C = 1+ 2 2 
I + 4, I being the orbital angular momentum contributed 
by the interacting pair, /j^ is a constant representing 
the overlap of the PH for interacting partition with the 
full set of all partitions, which can be found in Ref. 0. 
The correlated potential matrix element Vkk 1 (r) is given 
by 



Vkk 



l + z 



(r) = (hfhf,)^ J x {Pk P {z)v[t 

PtfWv ( r V^T ) Wt(z)}dz. (10) 



Here and Wi(z) are respectively the norm and 

weight function [ll| of the Jacobi polynomial P^?{z). 
Note that the inclusion of the short-range correlation 
function, r\{tij) makes the PH basis non-orthogonal. 
Numerical solution of Eq. © shows that r](rij) dif- 
fers from a constant value only in a small interval of 
small Tij values. Hence the dependence of the overlap 

<^aS+l("w ) )l^2^(^ )»7(» , H) > on the hypcrradius 
r is quite small. Disregarding derivatives of this overlap 
with respect to the hyperradius, we approximately get 
Eq. ©, with V K K'(r) given by Eq. dTJJJ) . The effect 
of the overlap being different from unity is represented 
by the asymptotic constant C of ^(r^ ). The emerging 
physical picture is: the effective interaction between 
pairs of atoms at very low energy becomes V(rij)r](rij). 
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This is justified, since at very low kinetic energy, the 
atoms have a very large de Broglie wave length and do 
not approach each other close enough to " see" the actual 
interatomic interaction. In the limit of zero energy, the 
scattering cross section becomes 47r|a s | 2 and the effective 
interaction is governed by the s-wave scattering length 
a s , through the asymptotic form of r](rij). 

Introduction of the PH basis and inclusion of the 
short-range correlation function, referred to as the 
correlated potential harmonic (CPH) method, simplifies 
the many-body problem dramatically. A fairly fast 
computer code can solve Eq. © with upto 15000 
particles in the condensate. This technique has been 
tested against known results, both experimental ones 
and theoretical ones calculated by other authors, for 
repulsive as well as attractive condensates [f| 0, Gil • 



III. NUMERICAL PROCEDURE 



momentum I. The HAA has been tested for nuclear, 
atomic and molecular systems and shown to give better 
than 1% accuracy, even for the long-range Coulomb 
potential [l3[. In our case, the van der Waals potential 
has a shorter range and HAA is expected to be better. 
Moreover, in a BEC, the dominant confining harmonic 
oscillator potential is smooth and the corresponding 
hypcrradial equation is completely decoupled. Hence in 
a BEC, the HAA is expected to be far better. We tested 
this by solving the CDE, Eq. ([3]), with the interatomic 
potential for the ground state by the renormalized 
Numerov method flil [l5j . which is an exact numerical 
algorithm for solving a set of coupled differential equa- 
tions. The calculated exact ground state energies are 
(in o.u.) 948.6420,1174.5284,1277.8219,1460.3706 and 
1596.2611 respectively for N = 700,900,1000,1200 and 
1400. These compare very well with the corresponding 
HAA results: 948.0986,1173.9809,1277.2040,1459.6373 
and 1595.8163 respectively. The error is less than 0.06% 
in all cases. Thus we can safely use the HAA, which 
reduce the numerical complications to a great extent. 



A. Solution of coupled equations 



Although Eq. ([9]) can be solved by an exact numerical 
technique using the Numerov method, we adopt the 
hyperspherical adiabatic approximation (HAA) fl2l ]. 
which apart from simplifying the computations greatly, 
provides an effective potential in the hyperradial space, 
in which the condensate moves. This effective conden- 
sate potential provides a physical picture for the internal 
mechanism of the condensate. 

In the HAA, one assumes that the hyperangular mo- 
tion is much faster than the hyperradial motion, since 
the latter corresponds to the breathing mode. There- 
fore, one can solve the former adiabatically for a fixed 
value of r and obtain the solution as an effective poten- 
tial for the hypcrradial motion, as in Born-Oppcnhcimcr 
approximation. The hyperangular motion is solved by di- 
agonalizing the potential matrix Vkk'{t) together with 
the hyper-centrifugal potential [second term of Eq. (j9|)]. 
The lowest eigenvalue ujq (r) [corresponding eigen column 
vector being Xxo^lLis the effective potential for the 



hyperradial motion 
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Co(r) = 0. (11) 



B. Calculation of specific heat 

At a temperature T > 0, bosons are distributed in 
available energy levels E n i according to Bosc distribution 
function 



f(Enl) 



1 



(12) 



where /3 = 1/ksT and [i is the chemical potential. The 
latter is determined from the constraint that the total 
number of particles is N. Clearly, [i has a temperature 
dependence. The total number of bosons in the trap is 
fixed and at any temperature it can be written as 
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(13) 



n=0 1=0 



At a particular temperature T, [i is determined from the 
constraint Eq. (fT3")) . The total energy of the system at T 
is given by 



The third term is an ovcrbinding correction. Eq. (|11[) is 
solved by the Runga-Kutta method, subject to appropri- 
ate boundary conditions to get Er and the hyperradial 
wave function Cq(t). The many-body wave function can 
be constructed in terms of (o(r) and XKo( r ) [Hi- Total 
energy is obtained by adding the center of mass energy 
(1.5 o.u.) to Er. Energy levels, E n i, are characterized 
by the quantum numbers (n, I), where n represents the 
excitation quantum number for a given orbital angular 
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E(N,T) = Y,Y,(2l + l)f(E nl )E nl (14) 



n=0 1=0 



The specific heat at fixed particle number N is calculated 
using the relation 



C N (T) 



dE(N, T) 



dT 



N 



(15) 
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Using (|T2j) , (fl4)) . (|15p one can obtain the heat capacity 
as 
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n=0 l n =0 



(2l n + l)E n i n exp - /i)) 

(exp(/3(£„ z „ -M))-l) 2 
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where 

dpi 



E~ =0 E~ rl= o( 2i >" + l)(E mim - M)oxp(f)(E mim - M))(/(E mim ))^ 
TEJL E i ~ = (2ip + l)oxp(/3(B plp - K»(/(E pip )) 2 

For an ideal non-interacting bosonic gas containing N 
bosons in a three-dimensional isotropic harmonic well, 
the critical temperature T® is well defined • M remains 
equal to the energy of the single particle ground state 
for T < T® and start decreasing rapidly for T > Tj?. In 
the standard text book treatment [l[ the sums in eqs. 
(|13H14|) are replaced in the semi-classical approximation 
by integrals over energy, assuming a continuous energy 
spectrum. In a harmonic trap the energy spectrum is 
discrete and this assumption is not valid, particularly for 
small N at low energies. A correct treatment Q shows 
that \x decreases slowly from its maximum value (equal 
to the ground state energy) as T increases from zero, 
the rate of decrease becoming suddenly rapid at some 
temperature close to the reference temperature T® corre- 
sponding to the same value of TV Q. Thus, in this case 
the critical temperature is not well defined. In the cor- 
rect treatment, the heat capacity also becomes a smooth 
function of T, attaining a maximum at a temperature, at 
which /i suddenly becomes a rapidly decreasing function 
of T. In the limit of large N, the behaviors of /Lt(T) and 
Cjv(T) curves approach those of the text book treatment. 
The transition temperature T c for a finite interacting sys- 
tem is defined as the temperature at which Cn{T) is a 
maximum [5| 

dC N (T) 



dT 







(18) 



For the numerical calculation for a chosen particle 
number N, the CPH equations are solved for a large 
number of energy levels - typically n running from to 
300 and / running from to 200, subject to an upper 
energy cutoff value, E UL (see below), so that E n % < E UL . 
Using previously calculated values of E n i, Eq. (fT3"|) is 
solved for fi, at a chosen temperature T, by a modified 
bisection method. Next Ejj^ is increased and the process 
repeated, until convergence in fi is achieved. Using this 
upper energy cutoff, Cjv(T) is determined using Eq. (Q~f 



IV. RESULTS AND DISCUSSIONS 

We consider the attractive condensate of 7 Li atoms in 
the trap used in the experiment at Rice University [l6j]. 



The harmonic trap used was axially symmetric with 
v x = v y — 163 Hz and v z = 111 Hz. For simplicity, 
we consider an isotropic trap with v = (y x v y v z )3. The 
experimental value of a s is —27.3 o.u. We use oscillator 

unit (o.u.) of length (y^j) and energy (frui). As 

mentioned earlier, we choose van der Waals (vdW) 
potential for the interatomic interaction, with known 
value of C 6 = 1.715 x 10~ 12 o.u. The value of r c is 
obtained by the procedure discussed following Eq. (Bl, 
so that calculated a s has the experimental value [3). 
Its numerical value is 5.338 x 10 -4 o.u. The calculated 
effective potential, u}q(t) is plotted as a function of r in 
Fig. [T]for N = 1300. In the r — > limit, uo(r) becomes 
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FIG. 1: (Color online) Calculated effective potential coo(r) 
against r in o.u. for the attractive 7 Li condensate with TV = 
1300. The narrow and deep well near origin is shown in the 
inset (note that different scales are used). 

strongly repulsive, due to the repulsive core of vdW po- 
tential and the hypercentrifugal repulsion of Eq. © . As 
r increases, there is a deep narrow well (DNW) arising 
from the strong interatomic attraction at small values of 
r. This attraction is proportional to the number of pairs 
and hence increases rapidly as N increases. For still 
larger r, the effects of the kinetic pressure (including the 
centrifugal repulsion), interatomic attraction and the 
harmonic confinement together produce a metastable 
region (MSR). An intermediate barrier (IB) appears 
between the DNW and MSR. The DNW is very deep 
and narrow, hence it is shown as an inset in Fig. [T] (note 
large changes in scale for both horizontal and vertical 
axes). As N increases, the DNW becomes deeper, IB 
shallower and the minimum of the MSR higher. At the 
critical value N cr , the maximum of IB and the minimum 
of MSR merge to form a point of inflexion and the MSR 
disappears. At this point, the condensate falls into the 
DNW, resulting in a collapse of the condensate and 
formation of clusters within the DNW. Our calculated 
value of N cr is 1430. In panels (a) - (f) of Fig. [21 we 
demonstrate how the MSR shrinks, with N approaching 
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FIG. 2: (Color online) Plot of effective potential uo(r) against 
r (both in appropriate o.u.) for the attractive 7 Li condensate 
with TV = 500, 1300, 1400, 1410, 1420 and 1426 in panels (a) 
- (f ) respectively, showing how the MSR shrinks in depth and 
width. Note that different scales have been used in different 
panels, to bring out the features of the MSR as TV — >• TV cr . 



N cr , for TV = 500,1300,1400,1410,1420, and 1426 re- 
spectively. From Fig. [5J one notices that both the depth 
and width of MSR decrease as TV increases towards N cr . 
Hence the number of bound energy levels supported by 
the MSR decreases rapidly with TV (see also Fig. [7J. 

However, the effective potentials shown in Fig. [TJ and 
Fig. [5] are obtained for I = 0. For higher I, the effective 
potential has a higher IB, arising from the /-dependent 
terms of the hyper-centrifugal repulsion [see Eq. (|5J)]. 
Thus the position of the MSR rises higher in energy as I 
increases, as can be seen in Fig. [3]for I = 0, 1, 2, 3, 4 for a 
condensate containing 1420 atoms. Hence particles with 
I > can attain higher energy levels. Inclusion of these 
levels will have a profound effect on the heat capacity. 
If such energy states were ignored (i.e., only the energy 
levels supported by the I = MSR considered), the 
heat capacity would reduce drastically and T c would 
increase indefinitely as TV — > N cr , since all the atoms 
would be forced into the few remaining energy levels 
available for internal excitation. In our calculation, we 
have retained all energy levels supported by a given I. A 
question arises as to whether the metastable condensate 
can have large I values. Intuition indicates that, with 
increase of temperature, the system can absorb energy 
only by increasing its rotational kinetic energy, thereby 
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FIG. 3: (Color online) Plot of effective potential lj 1 (r) against 
r (both in appropriate o.u.) for the attractive 7 Li condensate 
with TV = 1420 atoms, corresponding to I = 0, 1,2,3,4. The 
curves show how the IB increases as I increases. 



increasing the stability of the metastable system with 
enhanced centrifugal repulsion. Increase of kinetic 
energy due to faster linear motion alone would cause 
the system to fall in the DNW near the center of the 
condensate. Compared with the non-interacting or 
repulsive condensates, the attractive condensate has a 
clear distinction, viz. while the number of hyperradial 
excitations for a given I in the former is not limited, it is 
drastically limited in the latter. Thus for an attractive 
condensate, there are fewer energy states, in which the 
system can reside. This causes (C^) max to increase 
initially for TV <C N cr (when energy levels are not greatly 
restricted), but as TV — 5- TV cr , it starts decreasing, after 
attaining a maximum. Fig. [4] shows how Cn(T) depends 
on T, for selected values of N. It is seen that the 
transition temperature T c increases gradually with N, 
but {CN)max increases up to TV = 1300 and for larger 
N, it starts decreasing. When TV <C N C r, the nature is 
similar to that of a repulsive condensate [J] , since in this 
case, the number of available energy levels are still large 
enough (the top most energy level - including I =/= 
in the MSR has an energy much greater than fc^T), so 
that the top most levels are still practically unoccupied 
and there is scope for further internal excitation as T in- 
creases. Consequently T c increases gradually with TV, as 
in the repulsive case. As TV approaches TV cr , the number 
of energy levels supported by the MSR decreases rapidly 
and there is less scope for absorbing energy internally as 
T increases. Hence (Cisi)max decreases and T c increases 
faster, as TV increases towards TV cr . At higher temper- 
atures, higher I states are excited, which push atoms 
further outwards, increasing the average interatomic 
separation. Consequently, the system behaves ultimately 
as a non-interacting Bose gas. Thus the asymptotic 
value of Cn{T) becomes 3Nks- In Fig. 0] we plot 
the dimcnsionlcss quantity Cn(T) / '(TVfcs) against T (in 
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FIG. 4: (Color online) Plot of heat capacity C N {T)/{Nk B ) 
(dimensionless) against T (in nK), for indicated number of 
7 Li atoms in the metastable condensate. 

nK) for 7 Li condensate with N = 500, 1000, 1300, 1350 
and 1400. The features discussed above are clearly 
visible. One notices that the behavior for N < 1300 
is similar to that of a repulsive condensate, but as 
N exceeds 1300, the curves become flatter near their 
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FIG. 5: (Color online) Plot of the peak heat capacity, 
(CN)max I (iVfes) (dimensionless), against the number N of 
bosons in the attractive 7 Li condensate. 

maxima and the peak value of CV(T) / '(Nks), namely, 
{CN)max/{Nk B ), decreases fairly rapidly, as N —> N cr . 
All the curves appear to converge to the Bose gas 
limit. But a closer scrutiny shows that the curves for 
TV = 1350 and 1400 show a slight downward trend. 
This is due to a limitation in the higher energy cut-off 
used in our calculation. In Fig. \E\ we plot calculated 
{CN)max/{NkB) as a function of N. It is seen that this 
quantity increases gradually up to N = 1300. Beyond 
this value, {C^) max j '(Nkg) decreases fairly rapidly as 
N — > N cr . A plot of transition temperature T c (in nK) as 
a function of N is shown in Fig. [5] Initially T c increases 



N 

FIG. 6: (Color online) Plot of transition temperature (in nK) 
versus N for the attractive 7 Li condensate. 

linearly for N < 1300. As discussed above, this behavior 
is expected for small iV, as in the case of a repulsive 
condensate. But for TV > 1300, T c increases rapidly. 
Both the decrease of (CN)max and faster increase of T c 
are due to reduction in the number of available energy 
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FIG. 7: (Color online) Plot of total number of available I = 
energy levels versus N, close to N cr for the attractive 7 Li 
condensate. 

levels as N approaches N cr . We demonstrate this in Fig. 
[7] for the I = energy levels as N increases from 1200 
to Ncr' The decrease in the number of available energy 
levels forces a larger fraction of the bosons to be in the 
ground state as T increases [see Eq. (|T3|) ]. This causes 
a decrease of (Cjy)mox an d an increase of T c , as N — > N cr . 

A possible scenario of the attractive condensate as its 
temperature is gradually increased is the following. The 
standard definition of critical number N cr referred to in 
the literature, corresponds to the I = condensate at 
zero temperature. As T is gradually raised, the system 
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absorbs energy by occupying higher available energy lev- 
els upto the top of the MSR. However, there is a finite 
life time of atoms in higher energy levels due to tunnel- 
ing through the IB into the DNW. Thus there will be a 
decrease in the number of atoms in the MSR. The rate of 
loss of atoms will increase with the energy of the level, as 
also with N (increase of N will lower the IB). The usual 
definition of heat capacity, Cn (T) , is the rate of change 
of internal energy with respect to T [Eq. (|15|) ]. for a 
fixed number N of atoms in the condensate. This defi- 
nition is unambiguous for a repulsive or non-interacting 
condensate, since in these cases there is no loss. However, 
since the loss is appreciable for an attractive condensate 
for N close to N cr , or if T is such that ksT is compa- 
rable with the highest excitation energy allowed by the 
IB, this definition demands that atoms be pumped into 
the condensate at the same rate as the loss rate from 
the condensate. Such an experimental procedure has not 
been adopted yet. However, for N <C N cr and T <C T c , 
the highest appreciably occupied levels will have negli- 
gible tunneling probability into the DNW. Under these 
conditions, the metastable attractive condensate is fairly 
long lived and the standard definition of Cn(T) is ac- 
ceptable. Hence standard experimental techniques can 
be adopted. Thus our results presented in Figs. 4-6 are 
experimentally verifiable in the small N, small T limit. 
We have presented, for theoretical completeness, results 
for N close to N cr and for T beyond T c . The question 
of how the system can absorb energy internally for such 
values of N and T was already discussed above. 



V. CONCLUSIONS 

In this work, we report a detailed calculation of the 
heat capacity Cn(T) of an attractive Bose-Einstein 
condensate containing N atoms of 7 Li. The correlated 
potential harmonics method, which is appropriate for 
the dilute BEC, has been used. The effective potential, 
in general, supports a large number of energy levels. 
At T = 0, the lowest energy level accommodates all 
the bosons. As temperature increases, particles are 
distributed in higher energy levels, according to Bose 
distribution formula. Thus the internal energy of the 
system increases. Heat capacity Cn for a fixed number 
N of particles in the condensate is defined as the 
temperature derivative of the total internal energy. For 
a repulsive condensate trapped by an ideal harmonic 
oscillator, the effective potential has no upper cut 
off. Hence the energy levels are not limited in energy. 
Consequently, total internal energy and Cn increase 
as temperature increases upto T c . For T > T c , the 
ground state occupation becomes suddenly microscopic 
(negligible) and the system behaves like a harmonically 
trapped Bose gas. Hence Cn decreases rapidly above 
T c , reaching its asymptotic value 3Nks- Thus Cn first 
increases, reaches a maximum value (CN)max and then 
decreases rapidly to its asymptotic value, as T increases 



from zero. 

For N < N cr bosons with mutual attraction, a 
metastable condensate is formed in the metastable 
region (MSR) of the effective potential. On the left of 
the MSR, an intermediate barrier (IB), followed by a 
deep narrow well and finally a strongly repulsive core 
appear, as one approaches the center of the condensate. 
For N <C N cr , the IB is very high and the minimum 
of the MSR is very low, so that the metastable well 
is sufficiently deep compared with thermal excitation 
energy ksT c at T c , and a large number of energy 
levels are supported. Hence for T < T c , even the most 
thermally excited particles do not feel the effect of the 
IB and Cn increases gradually, as in the repulsive case. 

With increase of temperature, the system with a fixed 
N absorbs energy internally by increasing the occupa- 
tion probability of higher energy levels supported by the 
metastable region of the effective potential. Atoms in 
energy levels close to the top of the intermediate barrier 
have appreciable probability to tunnel into the deep nar- 
row well, causing the condensate to loose atoms. But 
such levels are not occupied with any appreciable proba- 
bility if N <C N cr and T « T c . Hence such a condensate 
is quasi-stable and Cn(T) calculated for a fixed N is ap- 
propriate. When the rate of atom loss from the conden- 
sate is appreciable, standard definition of heat capacity 
at constant N requires feeding the attractive condensate 
with additional atoms at a rate such as to compensate for 
the loss rate. Although this is not the usual experimen- 
tal technique, we investigate such cases for a complete 
theoretical study. In such a situation, there are dras- 
tic changes. The peak value of CV(T) (the temperature 
at which this occurs is the transition temperature T c ) 
initially increases gradually with N , then after reaching 
a maximum, decreases fairly rapidly near N cr . On the 
other hand, for small N, T c increases almost linearly up 
to N ~ 1300. For larger N, the transition tempera- 
ture increases rapidly with N. We provide an explana- 
tion of this behavior, based on the microscopic mecha- 
nism of absorption of internal energy, as T increases. As 
N increases towards N cr , the depth and width of the 
metastable well decrease rapidly. As a result, the num- 
ber of energy levels supported by the metastable well 
decreases rapidly. This tends to increase T c , since fewer 
energy levels are available for absorption of internal en- 
ergy, and bosons are forced to be in lower energy levels as 
T increases. Rapid reduction of available energy levels as 
N — > N cr , causes quicker saturation of internal energy of 
the condensate. Consequently, the maximum of Cn(T) 
decreases rapidly as N — > N cr . 
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